home *** CD-ROM | disk | FTP | other *** search
/ Software Vault: The Diamond Collection / The Diamond Collection (Software Vault)(Digital Impact).ISO / cdr44 / newmat08.zip / TMTD.CPP < prev    next >
C/C++ Source or Header  |  1995-01-11  |  6KB  |  182 lines
  1.  
  2. //#define WANT_STREAM
  3.  
  4. #include "include.h"
  5. #include "newmatap.h"
  6.  
  7. void Print(const Matrix& X);
  8. void Print(const UpperTriangularMatrix& X);
  9. void Print(const DiagonalMatrix& X);
  10. void Print(const SymmetricMatrix& X);
  11. void Print(const LowerTriangularMatrix& X);
  12.  
  13. void Clean(Matrix&, Real);
  14.  
  15.  
  16.  
  17.  
  18. void trymatd()
  19. {
  20. //   cout << "\nThirteenth test of Matrix package\n";
  21.    Tracer et("Thirteenth test of Matrix package");
  22.    Exception::PrintTrace(TRUE);
  23.    Matrix X(5,20);
  24.    int i,j;
  25.    for (j=1;j<=20;j++) X(1,j) = j+1;
  26.    for (i=2;i<=5;i++) for (j=1;j<=20; j++) X(i,j) = (long)X(i-1,j) * j % 1001;
  27.    SymmetricMatrix S; S << X * X.t();
  28.    Matrix SM = X * X.t() - S;
  29.    Print(SM);
  30.    LowerTriangularMatrix L = Cholesky(S);
  31.    Matrix Diff = L*L.t()-S; Clean(Diff, 0.000000001);
  32.    Print(Diff);
  33.    {
  34.       Tracer et1("Stage 1");
  35.       LowerTriangularMatrix L1(5);
  36.         Matrix Xt = X.t(); Matrix Xt2 = Xt;
  37.         QRZT(X,L1);
  38.       Diff = L - L1; Clean(Diff,0.000000001); Print(Diff);
  39.         UpperTriangularMatrix Ut(5);
  40.       QRZ(Xt,Ut);
  41.       Diff = L - Ut.t(); Clean(Diff,0.000000001); Print(Diff);
  42.       Matrix Y(3,20);
  43.       for (j=1;j<=20;j++) Y(1,j) = 22-j;
  44.       for (i=2;i<=3;i++) for (j=1;j<=20; j++)
  45.          Y(i,j) = (long)Y(i-1,j) * j % 101;
  46.       Matrix Yt = Y.t(); Matrix M,Mt; Matrix Y2=Y;
  47.       QRZT(X,Y,M); QRZ(Xt,Yt,Mt);
  48.       Diff = Xt - X.t(); Clean(Diff,0.000000001); Print(Diff);
  49.       Diff = Yt - Y.t(); Clean(Diff,0.000000001); Print(Diff);
  50.       Diff = Mt - M.t(); Clean(Diff,0.000000001); Print(Diff);
  51.         Diff = Y2 * Xt2 * S.i() - M * L.i();
  52.       Clean(Diff,0.000000001); Print(Diff);
  53.    }
  54.  
  55.    ColumnVector C1(5);
  56.    {
  57.       Tracer et1("Stage 2");
  58.       X.ReDimension(5,5);
  59.       for (j=1;j<=5;j++) X(1,j) = j+1;
  60.       for (i=2;i<=5;i++) for (j=1;j<=5; j++)
  61.          X(i,j) = (long)X(i-1,j) * j % 1001;
  62.       for (i=1;i<=5;i++) C1(i) = i*i;
  63.       CroutMatrix A = X;
  64.       ColumnVector C2 = A.i() * C1; C1 = X.i()  * C1;
  65.       X = C1 - C2; Clean(X,0.000000001); Print(X); 
  66.    }
  67.  
  68.    {
  69.       Tracer et1("Stage 3");
  70.       X.ReDimension(7,7);
  71.       for (j=1;j<=7;j++) X(1,j) = j+1;
  72.       for (i=2;i<=7;i++) for (j=1;j<=7; j++)
  73.          X(i,j) = (long)X(i-1,j) * j % 1001;
  74.       C1.ReDimension(7);
  75.       for (i=1;i<=7;i++) C1(i) = i*i;
  76.       RowVector R1 = C1.t();
  77.       Diff = R1 * X.i() - ( X.t().i() * R1.t() ).t(); Clean(Diff,0.000000001);
  78.       Print(Diff);
  79.    }
  80.  
  81.    {
  82.       Tracer et1("Stage 4");
  83.       X.ReDimension(5,5);
  84.       for (j=1;j<=5;j++) X(1,j) = j+1;
  85.       for (i=2;i<=5;i++) for (j=1;j<=5; j++)
  86.          X(i,j) = (long)X(i-1,j) * j % 1001;
  87.       C1.ReDimension(5);
  88.       for (i=1;i<=5;i++) C1(i) = i*i;
  89.       CroutMatrix A1 = X*X;
  90.       ColumnVector C2 = A1.i() * C1; C1 = X.i()  * C1; C1 = X.i()  * C1;
  91.       X = C1 - C2; Clean(X,0.000000001); Print(X);
  92.    }
  93.  
  94.  
  95.    {
  96.       Tracer et1("Stage 5");
  97.       int n = 40;
  98.       SymmetricBandMatrix B(n,2); B = 0.0;
  99.       for (i=1; i<=n; i++)
  100.       {
  101.          B(i,i) = 6;
  102.          if (i<=n-1) B(i,i+1) = -4;
  103.          if (i<=n-2) B(i,i+2) = 1;
  104.       }
  105.       B(1,1) = 5; B(n,n) = 5;
  106.       SymmetricMatrix A = B;
  107.       ColumnVector X(n);
  108.       X(1) = 429;
  109.       for (i=2;i<=n;i++) X(i) = (long)X(i-1) * 31 % 1001;
  110.       X = X / 100000L;
  111.       // the matrix B is rather ill-conditioned so the difficulty is getting
  112.       // good agreement (we have chosen X very small) may not be surprising;
  113.       // maximum element size in B.i() is around 1400
  114.       ColumnVector Y1 = A.i() * X;
  115.       LowerTriangularMatrix C1 = Cholesky(A);
  116.       ColumnVector Y2 = C1.t().i() * (C1.i() * X) - Y1;
  117.       Clean(Y2, 0.000000001); Print(Y2);
  118.       UpperTriangularMatrix CU = C1.t().i();
  119.       LowerTriangularMatrix CL = C1.i();
  120.       Y2 = CU * (CL * X) - Y1;
  121.       Clean(Y2, 0.000000001); Print(Y2);
  122.       Y2 = B.i() * X - Y1; Clean(Y2, 0.000000001); Print(Y2);
  123.  
  124.       LowerBandMatrix C2 = Cholesky(B);
  125.       Matrix M = C2 - C1; Clean(M, 0.000000001); Print(M);
  126.       ColumnVector Y3 = C2.t().i() * (C2.i() * X) - Y1;
  127.       Clean(Y3, 0.000000001); Print(Y3);
  128.       CU = C1.t().i();
  129.       CL = C1.i();
  130.       Y3 = CU * (CL * X) - Y1;
  131.       Clean(Y3, 0.000000001); Print(Y3);
  132.  
  133.       Y3 = B.i() * X - Y1; Clean(Y3, 0.000000001); Print(Y3);
  134.  
  135.       SymmetricMatrix AI = A.i();
  136.       Y2 = AI*X - Y1; Clean(Y2, 0.000000001); Print(Y2);
  137.       SymmetricMatrix BI = B.i();
  138.       BandMatrix C = B; Matrix CI = C.i();
  139.       M = A.i() - CI; Clean(M, 0.000000001); Print(M);
  140.       M = B.i() - CI; Clean(M, 0.000000001); Print(M);
  141.       M = AI-BI; Clean(M, 0.000000001); Print(M);
  142.       M = AI-CI; Clean(M, 0.000000001); Print(M);
  143.  
  144.       M = A; AI << M; M = AI-A; Clean(M, 0.000000001); Print(M);
  145.       C = B; BI << C; M = BI-B; Clean(M, 0.000000001); Print(M);
  146.  
  147.  
  148.    }
  149.  
  150.    {
  151.       Tracer et1("Stage 5");
  152.       SymmetricMatrix A(4), B(4);
  153.       A << 5
  154.         << 1 << 4
  155.         << 2 << 1 << 6
  156.         << 1 << 0 << 1 << 7;
  157.       B << 8
  158.         << 1 << 5
  159.         << 1 << 0 << 9
  160.         << 2 << 1 << 0 << 6;
  161.       LowerTriangularMatrix AB = Cholesky(A) * Cholesky(B);
  162.       Matrix M = Cholesky(A) * B * Cholesky(A).t() - AB*AB.t();
  163.       Clean(M, 0.000000001); Print(M);
  164.       M = A * Cholesky(B); M = M * M.t() - A * B * A;
  165.       Clean(M, 0.000000001); Print(M);
  166.    }
  167.    {
  168.       Tracer et1("Stage 6");
  169.       int N=49;
  170.       int i;
  171.       SymmetricBandMatrix S(N,1);
  172.       Matrix B(N,N+1); B=0;
  173.       for (i=1;i<=N;i++) { S(i,i)=1; B(i,i)=1; B(i,i+1)=-1; }
  174.       for (i=1;i<N; i++) S(i,i+1)=-.5;
  175.       DiagonalMatrix D(N+1); D = 1;
  176.       B = B.t()*S.i()*B - (D-1.0/(N+1))*2.0;
  177.       Clean(B, 0.000000001); Print(B);
  178.    }
  179.  
  180. //   cout << "\nEnd of Thirteenth test\n";
  181. }
  182.